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Abstract. We investigate how to coherently define entropy production for a process of 
transient relaxation in the Quantum Brownian Motion model for harmonic potential. We 
compare a form, called "Poised" (P), which after non-Markovian transients corresponds to 
a definition of heat as the change in the system Hamiltonian of mean force, with a recent 
proposal by Esposito et al (ELB) based on a definition of heat as the energy change in the 
bath. Both expressions yield a positive-defined entropy production and coincide for vanishing 
coupling strength, but their difference is proved to be always positive (after non-Markovian 
transients disappear) and to grow as the coupling strength increases. In the classical over- 
damped limit the "Poised" entropy production converges to the entropy production used in 
stochastic thermodynamics. We also investigate the effects of the system size, and of the 
ensuing Poincare recurrences, and how the classical limit is approached. We close by discussing 
the strong-coupling limit, in which the ideal canonical equilibrium of the bath is violated. 



1. Introduction 

The theory of stochastic thermodynamics provides a consistent description of nonequihbrium 
thermodynamics for classical systems weakly coupled to their environments and described 
by Markovian dynamics [HE]- In recent years fundamental characteristics of classical 
thermodynamics have been put under scrutiny in the quantum realm where a proper formulation 
of nonequihbrium thermodynamics seems a much harder task. In particular when considering 
low temperatures and non- vanishing couplings various difficulties arise, some of which are already 
present at equilibrium. An ubiquitous exactly solvable model to address these questions is the 
Quantum Brownian Motion (QBM) model. It consists of a system with Hamiltonian Hs (often an 
harmonic oscillator) bi-linearly coupled via a term denoted Hj to a bath of harmonic oscillators 
with Hamiltonian Hb- The total Hamiltonian is thus of the form H = Hs + Hb + Hj. In this 
model, when the total system is initially in canonical equilibrium p'^'^ = e~^^ /Z, the Clausius 
formulation of the second law seems to be violated for a quasi-static change of the mass or of 
the frequency of the central oscillator [3l[3]. In this case, the heat flow is defined as the change 
in the averaged central system Hamiltonian due to the bath and is found to be larger than the 
temperature times the change in the system entropy, defined as the von Neumann entropy of 
the central system. To remain consistent with these definitions, work and free energy are also 
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defined in terms of the system Hamiltonian and as a result tlie Thompson formulation of the 
second law is also violated E] • More work can be extracted from the system than the change 
in its free energy. One also intriguingly finds that, when the total system is in p^'^, the behavior 
of the heat capacity of the system at low temperature depends on its definition [71 [8]. Also, 
the von Neumann entropy of the central oscillator does not vanish at zero temperature while 
the equilibrium entropy of the total system (which coincides with the von Neumann entropy of 
the total system) does. These phenomena can be ascribed to quantum correlations between the 
central system and the bath |9]. Various attempts have been made in order to overcome these 
difficulties. Some of them incorporate in Thompson's formulation of the second law the work 
contribution required to couple the system to the bath at zero [5] or arbitrary temperatures 
|10| [TT] . Others introduce different notions of effective temperature [6l [12] . Ultimately, many of 
the difficulties are related to the fact that the equilibrium density matrix of the central oscillator 
is not the familiar canonical distribution = e~^^^ /Zg defined in terms of the central system 
Hamiltonian as is often the case in weak-coupling theories. 

In this paper, we want to investigate the slightly different problem of transient relaxation 
to equilibrium in the QBM model. This means that we initially place the central oscillator 
in a nonequilibrium state and put it in contact with its bath at equilibrium. Due to the 
interaction, the two parts of the system will exchange energy and if the bath is sufficiently large, 
the central oscillator will asymptotically reach an equilibrium state. For such a process we would 
like to identify a meaningful notion of entropy production. In stochastic thermodynamics the 
nonequilibrium version of the second law states that for such a relaxation process the entropy 
production is equal to the change in system entropy (which is identified with the Shannon 
entropy of the system) minus the heat exchanged with the bath divided by the bath temperature. 
Furthermore the entropy production can be proved to be an always positive quantity which only 
vanishes at equilibrium. For quantum systems the Shannon entropy is replaced by the von 
Neumann entropy of the system 5 = — trs ps In ps ■ An attempt to define entropy production for 
the QBM model was proposed in [6] and relied on defining the heat exchanged with the bath 
in terms of the system Hamiltonian (by integrating Q = trgffgpg) but with non satisfactory 
results, since a negative entropy production rate was obtained for non-vanishing system-bath 
coupling strength. A more satisfactory definition of entropy production has been recently 
introduced by Esposito, Lindenberg and Van den Broeck |13] (denoted here by ELB) in the 
form AiS = AS — Q/T, where the heat is now defined as minus the energy change in the bath 
{Q = — ti HbPb)- This quantum entropy production is positive definite even for finite bath 
sizes, notwithstanding recurrences. This definition applies on the assumption that the central 
system and the bath have uncorrelated density matrices at the initial time. 

In this paper, we compare this definition of entropy production A^S with a new definition 
AiS^ , inspired by the one introduced in [14j in the context of Markovian master equations. 
The corresponding heat definition is expressed in term of the averaged change of an effective 
Hamiltonian which reduces to the system Hamiltonian Hs in the weak-coupling limit. We 
evaluate analytically both expressions of the entropy production for the QBM model and evaluate 
their difference. We consider only Gaussian initial conditions, both for the bath and for the 
central oscillator, what guarantees that the density matrix remains Gaussian at all times. We 
find that A-^S is positive definite but can present oscillations while Ai5^ is positive definite and 
has a positive time derivative only in the Markovian high-temperature or weak-coupling limits. 
The difference between the two definitions sensibly depends on the coupling. We also study 
the behavior of the entropy production for finite-size thermal baths, where Poincare recurrences 
characterize the time evolution of the system. The convergence towards a continuous relaxing 
behavior is studied as a function of the system size. It turns out that a Lorentzian, rather than 
uniform, sampling of oscillation frequencies of the bath guarantees a better convergence. Finally 
the evolution the von Neumann entropy of the bath is studied. It is found that, for fixed initial 



conditions, its asymptotic value does not depend on the coupling in the classical limit while it 
does in the quantum regime. However, in both cases the Kullback-Leibler divergence between 
the density matrix of the bath at time t and at canonical equilibrium depends considerably on 
the coupling, making the usual approximation of the ideal bath problematic. 



Outline 

In section [2] the Quantum Brownian Motion model is introduced and solved. Our protocol is 
defined in subsection 2.1 and the evolution of the system is described in section 2.2 via its Wigner 
quasi-distribution function. The approach to the thermodynamic limit is described in section |3j 
The different definitions of the entropy production are spelled out in section [4j while their 
explicit expressions for our model are reported in section [5] Section [6] is devoted to a study 
of the model with a finite-sized bath, where Poincare recurrencies characterize its behavior. 
A discussion of the bath entropy, correlation entropy and of the distance of the bath density 
operator from its canonical form is reported in section[7j In section|8]we conclude and summarize 
our results. A few technical details are relegated in several appendices. 



2. The model 

The QBM Hamiltonian represents an harmonic oscillator bi-linearly coupled with coupling 
constants Ci to a bath of harmonic oscillators: 

H = Hs + Hb + Hi; (1) 
1 1 ^ ^ 

1=1 i=l 

We have put all masses equal to one for simplicity. 

The equations of motion in the Heisenberg picture read 

Q,it) = ^[H,Q^{t)] , P,{t) = ^[H,P^{t)]. (2) 

The Greek indices fx and u include by convention also the central oscillator and the terms 
associated to the bath, while the Latin ones only run on the bath degrees of freedom. For 
convenience, we will use the shorthand notation Q = Qo and P = Pq. The solution of the 
equations of motion reads [l5l [T6] 

N 

It is obtained by first finding a matrix transformation into new conjugate operators {Q'f^,Pl^} 
which diagonalize the Hamiltonian into a set of A^ + 1 normal harmonic oscillators, by then 
writing the Heisenberg solutions in that basis, and by finally transforming back to the original 
operators. The functions Anu{t) can be expressed in terms of the function 

N 2 

g(z) = z'-ul-Y,^^, (4) 

1=1 ' 

whose zeros z,^, v = 0, . . . , N , are the normal frequencies of the harmonic oscillators in the new 
basis. We have in fact ^ 

Aij,u{t) = Xfj,pX,yp " ^ ^ \ (5) 

p=0 



where the elements X^^ of the transformation matrix are given by 



Xi,, - 



"l dg{z) 




2z dz 


Z = Zi, 



-1/2 



e,; 



z^ = 0, 



(6) 
(7) 



From Q one gets also a condition to be fulfilled in order to obtain a positive-definite 
Hamiltonian and non diverging solutions: 



02 - , ,2 
J Zq — Wq 



$7o actually defines a normalized frequency of the central oscillator, as it emerges in the 
following Quantum Langevin Equation (QLE) (10). The role of the term A(t) = Aqo is that of 
a retarded propagator (see y^)- To see it, we first exploit the explicit solution to write down 
the equations of motion for the position operators in the bath as integro-differential equations 
involving the position operator of the central oscillator: 



Qiit) = Qi(0) Q.o's.iuJit) + ^''^^^ sin(wjt) - — 



As sin[a;j(t — s)\Q{^s). 



(9) 



Then the central oscillator satisfies the following QLE which we express in a matrix 
representation: 

z(t)^H^z{i) = -r](t)-F(t). (10) 
In this expression, * represents the time convolution, (i) is given by 







-m 

^ib(t) K{t) 



and we have defined z{t) with z'^{t) = {Q{t), P{t)) and the noise r]'^{t) = {0,ri{t)) with 
components 



TV 
i 

The damping kernel K{t) is defined by 

K{t)- 

where 7(w) is the coupling strength 



Qi (0) cos ujit H — sin Uit 



7(w) 



COSOjt, 



(12) 



(13) 



(14) 



and F(t) = (0, -fC(t)(5(0)) is the forcing term, which is responsible in the continuum limit for a 
fast slip of the initial conditions (see section [3|). 

The solution of the differential equation (10) can be easily obtained by taking the Laplace 
transform and then transforming back. One obtains 

z(t) = *(t)z(0)-(**T7)(t). (15) 



Both terms contain the matrix propagator ^(t) 



A{t) A{t) 
A{t) A{t) 



(16) 



which depends on the damping kernel K(t) via the propagator A(t). The Laplace transform of 
A{t) is given by 

^^""^ ~ s2 + sK{s) + ' ^^'^^ 



In Appendix A we show that the Ullersma solution ([s]) and the Fleming one ( 9|15 ) are 



equivalent. 

Having determined the time evolution of the Heisenberg momenta and positions as functions 
of the same operators at time t = 0, all the moments of these quantities at time t can now be 
evaluated as functions of the moments at t = and of the ^^i/'s. 

2.1. Initial conditions 

We assume to prepare the central oscillator and the bath separately, so that no correlation is 
initially present between them, and to instantaneously switch on the interaction at t = 0. We 
also assume that the bath is initially at canonical equilibrium. The density matrix of the total 
system is therefore of the form 



/o(0) = Ps(0)C5pb(0) , pb(0) = /9^^= , ZB = trBe-^^B^ (18) 

where /3 = {kBT)~^ is the Boltzmann factor and ps(0) and /Ob(0) are respectively the central 
oscillator and the bath reduced density matrix. From now on we set A;b = 1- 

An equivalent description of the system can be obtained via the Wigner quasi-probability 
distribution (often simply called "Wigner"), a function of the phase-space variables {q,p) = 
{qo,po, . . . ,qiy,PN), defined in term of the total density matrix p by 

W{q,p) = ^ rid^. (.0 - f , . . . , - f I P |.0 + |, . . . , + f ) . (19) 

The reduced Wigner corresponding to the central oscillator and the bath reduced density matrix 
can be defined in a similar way and can be obtained from the total system Wigner by integrating 
out the appropriate degrees of freedom. 

Using a matrix formalism with vectors z""" = {p, q), = (Q, P) and = {kg, kp), a generic 
single-particle Gaussian Wigner and its Fourier transform read 



N 1 r {z — {z)Y aiz — (z)) ^ , , 

W{k) = fi-x.Y>{-\k^c^k--\k^ {z)}. (21) 

The first moments and the symmetric covariance matrix respectively read {z)^ = ((P) , (Q)), 
{z)^ = {{Q),{P)) and <Tjj = {{zi, Zj}) /2 — {zi) (zj) with i,j = 1,2, where we denote by 
{...,...} the anticommutator and by (. . .) the average over a Gaussian density matrix p. We 
also define 

A = {ay^ - C'J^ ={det ct)K (22) 



where we indicate = an, = (T22 and Cqp = <ti2. 

We are going to restrict ourself to initial Wigner of the central oscillator which have a generic 
Gaussian form. It is then parametrized by the moments (Q(0)), (-P(O)) <|(5^(0)), (P^(0)) and 
Cqp{0) = ({(5(0) - (Q(0)) ,P(0) - (P(0))}) /2. Since the initial density matrix of the bath is 
a product of exponentials of quadratic Hamiltonians its corresponding Wigner is a product of 
Gaussian states which are parametrized, for i = 1, . . . , A^, by the moments 



(QiiO)) = {Pi{0)) = {{QM,P^m) = 



(23) 



where 



^, hw , hoj 

E{u;,T) = —coth—. 



As a result the initial total Wigner is also Gaussian. 



(24) 



2.2. Evolution 

At time t > the total density matrix operator will evolve in a unitary way: 

p{t) = u{t) {ps{o)^PBmuHt), 



-iHt 



(25) 
(26) 



and the central oscillator and the bath will be correlated. The corresponding total Wigner 
satisfies the Liouville-like evolution equation [161 EH] 



d_ 
di 



W{q,p,t) = [H,W] 



PB ' 



(27) 



where -ff ([T]) is now considered as a function of the phase-space variables {q,p) = 
{qo,Po, • • • ) Qn,Pn) and where [...,.. .]pg are the Poisson brackets. Again, by the linearity of the 
dependence of the solution for {q,p) on the initial conditions, an initial Gaussian distribution 
remains Gaussian at later times. This means that the Wigner is a real Gaussian, positive definite 
at all times, and fully characterized by its first and second moments. This also applies to the 
bath and will be useful to evaluate its entropy as we are going to see in section [7j 

It was shown in refs. |161I17| that the reduced Wigner satisfies the following partial differential 
equation: 

^Ws{z,t)= [vl-n{t)-z + Vl-D{t)-V,] Ws{z,t), (28) 
where the pseudo-Hamiltonian 'H{t) and diffusion D[t) matrices are respectively given by 



m 







-1 



nlit) 2T{t) 



-*(t)*-^(t), 



1 

'2^ IP 



1 r 
2 



n{t)aT{t) + aT{t)n^{t) + &T{t) 



(29) 



(30) 



where ^(t) is the matrix propagator (16) and the thermal covariance matrix (TT{t) is defined by 
^2 



a 



C, 



C, 



gp,T 



qp,T <^p,T 



dr [ dr' ^{t-r) 
Jo 




z^(r - r') 



(31) 



The noise kernel is defined as 



so that we have exphcitly 



dw 



7(0;) 



E{u}, T) cosujt, 



(32) 



CTT{t) 



duj 



E{u,T) 



1 d_ 

2 dt 



/o*(it'A(t')e- 



1 d_ 

2 dt 



/o dt' A{t')e'^'' 



(33) 



'H^t) and D{t) depend on the couphng strength via the damping and noise kernels. The 
off-diagonal elements of the diffusion matrix are called the anomalous dijfusion coefficients. 



The solution of equation (28) can be found by a Fourier transformation, via the method of 
characteristics [T7] : 

Ws{k, t) = Wsi^'^{t)k, 0) e-lfc^-^TWfc^ (34) 

which appears as a product of a function depending on the Wigner at time zero Wsik, 0), times 
a Gaussian one containing the thermal covariance. It clearly assumes Gaussian form in our 



hypotheses where the initial Ws{k, 0) is Gaussian (21 ). The dynamics of the central oscillator is 



then fully described by the first and second moments of the position and momentum operators: 



{zit)) = *(t)zo; 

a{t) = *(t)cro*"^(t) +crT(t). 



(35) 
(36) 



The general covariance matrix <T(t) corresponds to the covariance matrix in (20) if the averages 



(. . .) are evaluated with the total density operator at time t (25). It appears as the sum of the 



contribution of the evolution of the initial conditions and of the thermal covariance. The same 
expressions can be found by taking the average over initial conditions of the operators in the 
Heisenberg form ([sL expressing the Ajo(t)'s as functions of A{t) ( |A.l ) and by then using the 
coupling strength (14) [16]. The elements of the correlation matrix are reported in more detail 



m 



Appendix B 



At finite sizes one expects oscillatory behavior both for the dissipation and the diffusion 
coefficients. As in Ref. jT5], in the following we are going to assume an Ohmic form with a 
large cut-off for the couplings e? in (14). This choice enables us to obtain time-independent 



dissipation coefficients in the continuum frequency limit, while the diffusion ones only become 
time-independent in certain limits such as the high temperature limit. In general, however, this 
would not be the case: by assuming for example a sub-Ohmic coupling with a slower decay 
for larger frequencies, one would have time-dependent and nonlocal dissipation and diffusion 
coefficients at all times, even in the high-temperature limit |17j . 

3. Thermodynamic limit 

The thermodynamic limit of an infinite number of bath oscillators is obtained by substituting 
a continuous function 7(0;) to the discrete coupling strength (14). We choose the Drude-like 
Ullersma coupling strength [16\ [T5] 



7(0;) 



2 Ka^uj"^ 



(37) 



The parameter k tunes the strength of the coupling, while the cut-off a, which is introduced in 
order to eliminate ultra-violet divergences, can be associated to the bath memory time. In fact 



the damping kernel (13) with this coupling strength is given by 

K{t) = na e" 



-at 



(38) 



and thus decays over times of order a~^. 

The general form of the propagator A{t) foUowing from the strength (37) can be found in 
|16| . It is characterized by three time scales: il, T and A, deriving from the poles A, TziziO, of the 
Laplace transform of the propagator (17). They are obtained by solving the following equations 
in which the bare central oscillator frequency lvq, the coupling k and the cut-off a appear: 



A + 2r 



a 



' 



A 



a 



UJn 



(39) 



In this situation, the renormalized frequency Qq (jSl) of the central oscillator is simply given by 



"0 



na. 



(40) 



The propagator A{t) describes a noisy damped oscillator, where Vt is the characteristic frequency 
and A and F characterize the damping rates. When the time scale 1/A is much shorter than 



1/r and l/ri, the damping kernel K{t) ([TSj) becomes delta-like and the evolution equation (10) 
becomes local in time. This is obtained by taking the large cut-off limit, defined by 



We can then approximate ( 39 ) by 



A ~ a , r ~ k/2 
and the propagator A{t) assumes the form 

A\oc{t) = 7^sin(rit)e 



1 



-Tt 



(41) 



(42) 



(43) 



which is typical of a damped Ornstein-Uhlenbeck process. One can observe in fact that 
the pseudo-Hamiltonian H{t), at 0(l/a) and for t ^ 1/a, becomes equal to the the time- 
independent matrix 

" -1 " 



loc 



2T 



(44) 



Since f^o is real, the transition between real and imaginary Q. corresponds to the transition 
between the under-damped (F < Oq) and the over-damped regimes (F > Oq). The quantum 
and time-dependence features of our process are then contained only in the noise kernel (32), 

. / 1 



i>{t) = —Ka 
where we have defined 



cot ^e~"* + 



TT 



^ (aT/3)2 



h/3 

2^' 



We can thus define the following limits: 
The low-temperature limit: 

a > l/Tp] 

The high-temperature classical limit: 



(45) 



(46) 



(47) 



l/Tp > a; 



(48) 
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Figure 1. Transition between real (wliite) and imaginary (yellow) $7. In the large cut-ofF limit 



(41 ) it corresponds to the transition between under-damping F < r^o and over-damping F > JIq- 



For a/ujo < (on the left of the dotted line) one always has Jl^ > 0. 



The weak-coupling limit: 



F < n,i/Tp. 



(49) 



The noise kernel determines the thermal covariance matrix cr^ (31 ) and, via (30), the diffusion 
coefficients of the ME. Thus the quantum and time-dependence features will show up in these 
quantities. One can evaluate the covariance matrix by using the local propagator (43) inside the 



general expression (33), up to terms of 0(1 /a) 
one (see 



Using this propagator instead of the general 



\16\ eq.(7.10)]) does not affect either the covariance thermal matrix or the diffusion 
coefficients in the large cut-off limit, even for times t < 1/a, since only correction of 0(1/ a) 
arise [171. Since the coupling strength (37) is an even meromorphic function, the integrals 
appearing in (33) can be evaluated by a contour integration in the complex plane. 

Complete expression for the thermal correlation matrix in the large cut-off limit, which are 
exploited in the following for the calculation of the entropies, were derived in [16] and are reported 
in Appendix B Their quantum features are due to the presence of the function E{u},T), whose 



poles at u = ik 2t: /t^ with k any integer, give rise to the thermal transients, i.e., to terms which 
vanish on a time scale of order Tp. These terms are also responsible for the time-dependence of the 
diffusion coefficients of the master equation [TB]- In the high-temperature classical limit, where 
E{uj, T) approaches T, all the thermal transients vanish and the expressions of the covariance 
matrix simplify. 

One can deduce from equations (16), (34) and (43) that, since limt^oo ^loc (^) = 0, in the 
thermodynamic limit the system eventually loses all information on its initial conditions, and 
its distribution assumes the characteristic Gaussian form corresponding to the late-time thermal 
covariance matrix, as described in the next subsection. 

An important feature of this model is the presence of initial slips in the momentum average, 
in the non-thermal part of the averaged square momentum and of the correlation between Q 
and P. Using the local propagator (43) from t = 0, implies neglecting an initial evolution of the 
system during a short time of order 1/a, in which the central oscillator is subjected to an initial 



kick [TTldB]. One can easily observe, indeed, that 



i(0) = / iioc(t=0+) = -2r. (50) 

Our description will thus only be valid for t 2> 1/a. The effect of initial slips both on the 
moments and on the definitions of entropy production are discussed in [Appendix C 



Late-time covariance matrix 

In the thermodynamical limit it is possible to evaluate the long-time behavior of the diffusion 
coefficients and of the covariance thermal matrix. They are related by 

Dgp{^) = {P\^))-nl{Q\oo)); (51) 
Dpp{oo) = 2T {P\oo)) . (52) 

Thus the anomalous diffusion coefficients survive, since the right-hand side of the first equation 
of (51) does not vanish. This implies that equipartition does not hold in the general quantum 
case. 

Interestingly, as observed in [16], one obtains 

(g^(oo)> = (Q^>^^ , (P^(oo)> = (P^)^^, (53) 

namely that the stationary form of the central oscillator density matrix at t = cxd equals the 
traced canonical equilibrium one of the total system: 



ps(oo)=trBP^^ , , Z = tre-^^. (54) 



This does not mean of course that the total system equilibrates: p{oo) ^ p^^ [TO]. Furthermore, 
it has been shown in ^] that 

ps(oo) = p^g-^ ^ , Z^^ = tTse-^^s\ (55) 

where the equilibrium effective Hamiltonian H^'^ is given by 

= + Im^.u:!,Q^. (56) 

The effective frequency Ues and mass Mgs are respectively given by 



Expressions for \Q and \P /^^ can be found in [20] and [TB], and are reported in (B.15 



B.16). 



It is worth noticing that the traced canonical equilibrium density matrix (55) can be 
equivalently written in the form [3 [8] : 

-^M^ ' --^1^^^ ' "Z^ ^^^^ 

where an Hamiltonian of mean force has been introduced, which differs from the effective 
by the additive constant — /3(lnZg — InZg*^). Both partition functions have a well known 



analytical expression |20j. For the practical purpose of evaluating the Breuer entropy flow (70) 
we will use the effective Hamiltonian. However we emphasize that the use of the mean force 
Hamiltonian leads to exactly the same entropy production (since only the density matrix is 
involved) and heat flow (since only differences in energies are considered). 

The identities ( 53|54|55 ) do not generally hold in open quantum systems. They are however 
an important feature of our bilinear model and hold independently of the choice of the continuous 
limit strength. 



4. Entropy production 
4.I. The 'Poised' definition 

Let us consider a central system S whose density matrix ps{t) evolves according to the equation 

ps{t) = V{t)ps{0) = tre t/(t)ps(0) (E) PB{0)U\t), (59) 

where the bath is not necessarily at thermal equilibrium and U (t) is the unitary evolution 
operator in the global system (26). Let us define the poised density matrix Pg(t) as the solution 
of 

Vit)plit) = pUt). (60) 



The explicit form of this density matrix in our model is derived in Appendix D and is given by 
(89). We can then introduce the following definition of the entropy production: 

Ai^p = [Dipsmiput)) - Dipsmm)] , m 

where D(. . . || . . .) is the Kullback-Leibler divergence, defined by 

D{p\\p') =tTplnp-tvplnp'. (62) 



One can prove that the expression ( |61| ) is positive definite as follows. From equation ( |59| ) we 
obtain 

Dipsimm) = D{v{t)psm\v{t)put)) 

= D (tiB U{t)psiO) (g) pB{0)U\t)\\ tre U{t)pl{t) pB{0)U\t)^ 

< D{psm\pm), (63) 



where we have used the property of the Kullback-Leibler divergence 



D{pi\\p2)>D{tvBPi\\trBp2)- (64) 

The evolution of the system density matrix can be formally written as 

ps{t) = C{t)ps{t) , C{t) = V{t)V-\t), (65) 

where the operator C{t) in general depends on time. In this case, the evolution operator V{t) 
may be written as V{t) = TexplJ^ (ir£(r)} in which T indicates the time-ordering operator. 
However, in the Markovian case, when C{t) can be assumed time-independent, the poised density 
matrix pg coincides with the stationary one pg*, defined by Cp'^ = 0. In this case the poised 
density matrix is constant in time and, if it is unique, coincides with the stationary one, which 
in turn corresponds to the asymptotic form of the density operator 



(66) 



It follows from (63) and (66) that the time derivative of A^S^ is also positive: 



dt 



. j.^ D{psit + dt)\\pf)-D{psmpf) > Q 

dt-s-O dt ~ 



However, this result does not hold in general when C{t) depends on t. 

Because in our model pg = p^ ( 55 ) , following Breuer and Petruccione 
following definition of the entropy production: 



(67) 



we introduce the 



A,S^^ = Dipsmpl"") - Dips 



Pi') 



This expression has a straightforward physical interpretation, since it can be rewritten as 



AiS^' = AS- AeS^'i, 



(68) 



(69) 



i.e., as the difference between the entropy change AS, and the entropy flow AcS''^ identified by 
the variation of the averaged effective Hamiltonian of the central system: 



A,S 



Br 



I3A{HI^), (70) 
that in the classical high-temperature 



5.1 



where H^^ is defined in (56). We shall see in section 
limit and in the quantum weak-coupling limit this definition becomes equal to the one used in 
the usual stochastic thermodynamics setup. 



4-2. The ELB definition 

A definition of the entropy production which guarantees its positivity in all limits has been 
recently introduced by Esposito et al [T3j. It reads 



A,S = D{p{t)\\ps{t)0p'^) 



(71) 



This quantity is clearly positive definite and within our assumptions it expresses the second law 
in the form 

AiS = AS + PA{Hb) . (72) 
One has in fact, exploiting the conservation of the von Neumann entropy: 



D [pmpsit) p's') 



tr p{t) In p{t) - tr p{t) In ps (t) - tr p{t) In p"^ 

tr p^ In p^ + tr ps (0) In (0) - tr ps (t) In ps (t) - tr pB (t) In p"^ 

ASit) + [trp^^ {-PHb) - trpBit) (-PHb)] . (73) 



From ( 72 ) we observe that the entropy fiow 



A,S = -pA {Hb) = /3(A {Hs) + A {Hi) ), 



(74) 



is now identified as the change of the bath energy Hb times the inverse temperature f3, as usual 
in equilibrium thermodynamics. It is then proportional to the change of the central oscillator 
energy plus an interaction term which is discussed in section Sl The time derivative of ( 71 ) is 
not generally positive, as we will observe in 5.1.1 Since by (70) and (74) one has H'^^ — )• Hs 



and (Hi) — t- 0, we see that the two definitions coincide in the weak-coupling limit. We note that 
the positivity of the right-hand-side of ( 72 ) was also remarked in ref . |21] . 



4-3. Difference between the two definitions 
Since the following identity holds 



D{p{tWi pI^) = D{p{t)\\ps{t)r^pl^)+D{ps{t)\\pf), 



(75) 



we find, using (71 ) and (68), that the difference between the ELB definition Ai5 and the Breuer 
one AiS^'' is given by 



A,S - AiS^^^ = Dip{t)\\Ps ® Pb') - ^(Ps(O)llpf ). 



(76) 



Due to the unitary evolution of the total density operator p{t) ( 25|26[ ), one can recast the first 
term of the right hand side of ( 76 ) in the form 



By inserting this identity into (76) and noting that 



D{ps 



pf) = D{psiO)0p''^\\pf(E>pl^), 



we obtain 

AiS - AiS""' = -trps(O) (S) /9^^[lnC/t(t)p|t ® p^^f/(t) - In/^f p'^] 
Moreover due to the inequality 

D{ps{0)(^p''^\\U\t)p^i(S>pl^U{t)) 

> DitMpsio) ® pl^}\\tMuHt)pf PB^um) 

= D{psm\trB{UHt)pf pl^Uit)}), 



(77) 



(78) 



(79) 



it follows from (79) that 



AiS- AiS""' > trs PS (0) In V{t)p'i - In p'^ 



(80) 



(81) 



where we have introduced the evolution operator for the central system associated to the total 
adjoint dynamics, implicitly defined by 



V{t)p^^ =tTB{U\t)p^^ p'i.^Uit)} 



(82) 



In Appendix E we show that in our model, when the operator C{t) is time-independent the 
adjoint C{t) is also time-independent, and moreover that the latter identifies the same stationary 
density matrix, so that 

(83) 



Vit)pf 



Ps ■ 



Then the right-hand side of (81) vanishes and the ELB expression is strictly larger than the 
Breuer one: 

(84) 



AiS - AiS^' > 0. 



Over times sufficiently long for the dynamics to be Markovian the inequality (84) holds in any 
case, since equation (83) is always verified. 



5. Explicit forms of the entropy production 

We report here the expHcit forms of the entropy production, according to the P (61), the ELB 



(71) and the 'Breuer' expression obtained in (68). To evaluate them one needs to know the 



expressions of the entropy and the entropy flow. 

Since the central oscillator density matrix is gaussian at each time t, its von Neumann entropy 
entropy is given by [18] 



Sit) 



where we have defined 



trs ps(i)ln Psit) 

Ait) + 2) In Uit) + 2 



Ait) -I) In (Ait) -I 



Ait) = h-' ialit)alit) - Clit)) 



1/2 



(85) 



(86) 



which is a function of the correlation matrix at time t. One notices that Sit) is correctly defined 
if the uncertainty principle is satisfied. 

The 'Poised' entropy production can be written as 



A5- Ae^P; 

■i^ipsi^)-psimnp*sit). 



(87) 



The 'Poised' density matrix p*git) is gaussian with vanishing means of Q and P (as shown 
in 



Appendix D ) and is given by 
\np*sit) = - 



1 



In A*(t) 



A*(t) j^A*(i) + i 



where we have defined 



2A*(t) A* it) 



1/2 



1 ' 
2 



A*it) = h-' {al\t)a;\t)-Cl^\t)) 
A* it) = a;\t)Q' + afit)P' - C*,(t){Q, P}. 



The variances and correlation cr* p^(t) and C*p(t) are given in Appendix D 

{A*it))^-{A*{t)),^__A*it) + l 



A,S' 



2 A* it) 



In 



A* it) 



1 ' 
2 



(89) 

(90) 
(91) 

Then one obtains 
(92) 



where we have defined, for any operator O acting on the Hilbert space of S, (0)^ = tr psit)0. 
As for the Breuer entropy flow, AgS'^'' is straightforwardly given by the change in the effective 



energy Ha'^ (56). Then one only needs to know its second moments at time t, which in the 



finite case are obtained from the first moments and from the correlation matrix which appear 
in ( B.1|B.6 ), while in the continuum case one exploits the general expressions (33) (35) with the 
Ullersma coupling strength (37). The entropy fiow A^S is instead proportional to the change in 



the bath energy (^74|. To evaluate it, one needs rather to evaluate the average of the interaction 



energy term (Hi). By using the Ullersma strength (37) in the large cut-off limit, it turns out 
that for t » 1 /a 

(Hiit)) =Dgpit)-Ka(Q\t)). (93) 



This evaluation is reported in Appendix F Thus, by comparing definitions (74) and (70), we 



obtain the difference between the two entropy fiows is given by 



ASe - Ae^^"- = /3 [A (Hs) - A {H^'') - Ka {Q\t)) + D^pit)] 



(94) 



The difference of entropy production is the same with opposite sign. 

As aheady mentioned, the expressions we use in the continuum Umit for the three definitions 
of entropy only apply for t ^ 1/a, after the initial slip has taken place. Their contribution to 



entropy, which is reported in [Appendix C[ implies that the entropy flows and productions often 



do not start from 0, as one can observe in the following figures. 



5.1. The Markovian case 

Generally, in the limit of short-lived thermal transients, namely Fr^ <C 1, the generator of the 
dynamics C can be considered time-independent, since the diffusion coefficients (30) are close to 
their t = oo limit (51). As we have seen in section 4.1, in this case the entropy production 



definition (68) and its time derivative turn out to be consistently positive, as the system 
equilibrium density matrix does not depend on time, and the entropy flow is given by the 
average variation of the effective system Hamiltonian H^"^, as already observed in equation (70). 
This holds both in the high-temperature classical and weak-coupling limits [16]. 



5.1.1. Classical limit. In the high-temperature classical limit (48) all the quantum features 



of the system disappear and the anomalous diffusion coefficient vanishes, thus recovering 
equipartition (51), since (-P^)gq = T, and (Q^)gq = ^/^o- particular the equation satisfied 
by the Wigner has exactly the form of the Kramers equation for an oscillator in contact with a 
bath at temperature T [22l 



In this limit, since we have A(t) ^ 1, Vt, the system entropy 



( 85 ) assumes its classical form for a Gaussian distribution 



S{t) ~ l + lnA(t), 
Moreover the effective equilibrium energy is given by 



(95) 



(96) 



Thus the definition ( 70 ) of the entropy flow reduces to that of stochastic thermodynamics, which 



is defined as the average variation of the effective energy of the system, namely the classical 
one with the renormalized frequency VLq in place of ojq. This means that the corresponding 
definition of the entropy production coincides in this limit with the one introduced in the theory 
of stochastic thermodynamics for the Kramers equation [2lj . In the overdamped limit F ^ J7o 
the momentum equilibrates much faster than position and can thus be traced out. The entropy 
production assumes in this case the form proposed in the theory of stochastic thermodynamics 
for the overdamped Fokker-Planck equation [25j. As long as momentum has not yet fully 
equilibrated, the latter expression constitutes a lower bound to the former one since it results 
from a coarse graining procedure (see, e.g., [26]). Let us also note that by taking the weak 
coupling limit k — )• one gets in (96) the bare frequency wq, and that then the entropy fiow 



becomes exactly equal to the change in the central oscillator energy, divided by the temperature 
of the bath. 

We show in Figure [2] the high-temperature limit ( 48 ) of the difference between the different 
definitions of the entropy production with different coupling strengths. In the classical limit 



the anomalous diffusion term in (93) vanishes and the normal diffusion coefficient is time 



independent. This means that the Poised and Breuer expressions for entropy becom e eq ual: 
AiS^'' = AiS*^. Considering also the expression assumed by H'^ (96), the expression (94) for 
the difference of fiows simplifies to 



Ae5 - AeS^^-- = {{Q\t)) + (Q^(0)» 



(97) 




Figure 2. The different expressions for the entropy production (left) and the entropy flow 
(right), for different couplings k (top and center underdamped, bottom overdamped) in the 
classical regime. The parameters are: temperature T = 1000, a = 10, = Ij '7^(0) = 100, 
ct2(0) = 100, Cgp(O) = 10, (Q(0)) = 100, (P(0)) = 100. 



The same difference with opposite sign holds for the entropy production. In the classical limit 
the thermal part of {Q'^{t)) is proportional to 1/I3Qq. This means that the difference between 
the two expressions for the entropy production diverges, since for large k one has 

nl = uj^-Ka-^ 0. (98) 

This appears clearly in the figure, where the different expressions for the entropy production 
AiS and for the entropy flow A^S are shown for different coupling strengths k, both in the 
underdamped and the overdamped regime. 

The difference between the definitions is due to the fact that the expression AiS and the 
corresponding expression A^S of the entropy flow both diverge in the limit (98) as 

Ae5(oo) = 1 - 2^ - f K {Q'm + (^'(0)» • (99) 

However, the expression AjS^'' diverges only logarithmically like the von Neumann entropy: 

A5(oo) ~ I ln/3J]o| - S{0) (100) 



In fact the expression AgS^"^ does not diverge, since in the effective Hamiltonian (96) only the 
renormahzed frequency appears: Mefjcj^g — )• Qq. 

One notices in Figure [2] that both expressions of the entropy productions are positive, but 
that the ELB one, A5j, exhibits damped oscihations yielding a nonpositive time derivative. This 
can be directly seen from the fact that the time derivative of AiS'^'' is positive, due to the fact 
that the process (5.1) is time independent, and that the ELB one differs from it by a constant 
plus a term proportional to which is characterized by damped oscillations. 

We remark here that usually in literature the total Hamiltonian is renormahzed by a self 
interaction term, such that no positivity condition ([s]) has to be satisfied. In this case there 
would not be any divergence of which would be proportional to T/ujq but we could 

make the difference between the two definitions of entropy production arbitrarily large, as we 
could take k — )• oo |6j. 
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Figure 3. Entropy productions Ai5^, AiS^'^ and AiS at different values of the coupling k: 
0.01 (blue), 0.04 (green), 0.08 (red), 0.09 (cyan) with temperature T = 0.001 (left), and different 
temperatures T: 0.001 (blue), 0.01 (green), 0.1 (red), 1.0 (cyan), with k = 0.09 (right). Initial 
condition are fixed: ^^(O) = 1.0, ^^(O) = 1.0, Cgp(O) = 0, (g(0)) = 0, (P(0)) = 0. 



5.1.2. Weak-coupling limit Another case in which entropy flow is the same as the one 
defined in stochastic thermodynamics is the weak-coupling limit in the general quantum setting 
r <C ri, Tq'^ . Here things are more delicate than before, since the anomalous diffusion coefficient 



Dqp{t) is persistent to first order in the couphng F, just as the normal diffusion coefficient Dpp(t): 

Dqp{oo) = ^hVReliPil + Xtb) - + inrg)] + O (T'^) ; (101) 
vr 

Dpp{^) = E{n,T) + 0{r^), (102) 

where ipiz) is the digamma function. Anyway their contribution to (Q^)gq is different as Dpp{t) 
contributes to order one, while Dqp{t) to order T as seen by inverting (51). One gets then 
equipartition to first order in T: 

{Q')^^ = E{n,T)/n' + o{r) , (p^) = E{n,T) + o{r) , (103) 

where $7 can be approximated by to first order in T. This correspond to an equilibrium 
density matrix (55) corresponding to the equilibrium Hamiltonian 



^^s' = ^(^oQ' + i'') (104) 



which is the same as the central oscillator one ([T]). 



5.2. Low-temperature limit 



In the low-temperature limit (47) one expects that the Breuer entropy production expression 



(68), as well as its time derivative, can become negative. The Poised and ELB expressions (71) 
remain instead positive, while their time derivative can be negative. As we observe in Figure |3| 
for sufficiently low temperature and strong couplings, the expression AjS"^^ becomes negative, 
exhibiting an oscillatory behavior. At higher temperatures or weaker couplings the amplitude 
of the oscillations becomes smaller. Thus in these limits one obtains a positive definite entropy 
production, as well as a positive time derivative. We observe that the Poised and Breuer entropy 
production have the same asymptotic value, as expected, since /Og(oo) = p^. 



One notices that also in the low-temperature limit (47) the ELB expression can be orders 
of magnitude larger than the other two, due to the coupling term na {Q'^{t)) which appears in 
the entropy flow. This difference can be much larger respect to the classical case, due to the 
presence of the quantum terms contained in (|(5^(t)), which actually become more relevant than 
the classical one. 

6. Poincaire recurrences 

When the number N of bath oscillators is finite, the dynamics is characterized by a recurrent 
behavior, with a period identified by the Poincare recurrence time tp ~ 27r/ in.m{z^^i — z^) 
|15j . where the Zj^^s are the normal frequencies. We can interpret this recurrence as an almost 
periodic return to the initial decoupled state. Interestingly, while AiS remains positive by 
definition, one might have a negative AjS'^'', even in the classical case. When the size of the 
bath becomes larger, the recurrence time grows, and one expects that the entropy approaches 
its typical irreversible behavior, eventually relaxing to the equilibrium asymptotic value. 

In the present section we study this behavior in the two specific cases of uniform and 
Lorentzian frequency sampling, always assuming that the coupling strength converges to the 



Ullersma expression (37). Indeed, the density of states 6{uj — u}i) inside the coupling strength 
can be arbitrarily chosen. We evaluate the thermal covariance matrix components a'^rp[t), 
(Tpj,{t), and the equilibrium symmetrized autocorrelation function C(t), defined by 



C{t) = l{{Q{t),Qm) 



(105) 



We can also consider the Fourier transform of the correlation function C{t). Indeed, in the 
classical limit, the finite-size correlation function has the expression 



N ^2 

Civ(t) = T ^ cos(z,t). (106) 

We can thus represent the Fourier transform Cj^iuS) of C(t) by setting it equal to TXq^/ (z^ ^u), 
where Ajy = Zu — Zy-i^ and considering it as a function of = Zu- This quantity should approach, 
as A'^ —7- DO, the Fourier transform of C(t), which is given by 

(a;2-a;2)2(a2 + u;2) + «;2c,4 + 2«;a3(c^2_^2)- l-^^'^ 

We will see that the convergence to the large-size irreversible behavior is much slower for the 
uniform than for the Lorentzian sampling, and that, in the former case, the dynamics seems to 
remain characterized by underdamped oscillations even at large values of N. 

6.1. Sampling 

6.1.1. Uniform. The uniform sampling is obtained by considering N frequencies {i = 
1,2, ... , N) spaced by a constant A. The maximal frequency A will be denoted by Uc- The 
corresponding couplings are given by 



Q = jA-^— ^. 108 



Then the continuous-limit UUersma strength is obtained for N ^ oo, Uc ^ oo and A = Uc/N — >■ 
0. In this case the Poincare recurrence time is given by ~ 27r/A. 

6.1.2. Lorentzian. In order to obtain a faster convergence with longer Poincare recurrence 
times, and a better agreement with the continuum curve both in the under-damping and in the 
over-damping cases, one can adopt a Lorentzian sampling of frcqiiencies. Positive frequencies 
distributed with a Lorentzian density centered at a; = 0, with width oq are defined as 



U(, = ao tan 



.^+12 

with £ = 1. . .N and with the corresponding couplings 



(109) 



ei = JAe ^— ^, 110 

where = coi — w^-i, £ = 2,. . . ,N and Ai = oji. This sampling enables a high density of 
frequencies in the area around co = 0, then determining a long recurrence time. One can adjust 
the value of A^v in such a way as to have, for all values of N, 

e2 



J2;^,=na. (Ill) 



We shall refer to this case as the adjusted Lorentzian sampling. 



6.2. Results 

In Figure |4] we report the correlation function and the thermal part of second moments of 
the central oscillator in the classical continuum limit both for an under-damping and an over- 
damping set of parameters. These are compared with the results obtained in the finite case with 
= 600 bath particles, both with a uniform and Lorentzian sampling of the bath frequencies. 
The parameters ujc and oq are chosen so that the Ullersma's spectrum is sampled beyond the cut- 
off a, and the recurrence time is of the order of the characteristic relaxation time l/T. Finally 
curves obtained with the adjusted Lorentzian sampling are reported, where the parameter ag is 
chosen so that the recurrence time is much longer than l/F. 




Figure 4. Plot of C{t), 0"^ and Cp j-C^) in the continuous limit (dashed), compared with 
the corresponding curves obtained for N = 600 with a uniform frequency distribution with 
LOc = 30.0 (blue), a Lorentzian distribution of frequencies with ao = 20.0, (red) and an adjusted 
Lorentzian distribution with ao = 0.1 (green). They are obtained both for an under-damping 
set of parameters (left column) and an over-damping one (right column) in the classical case 
T = 1000. Insets are magnifications of the finite-size curves with uniform sampling. 

In the under-damping case, for a finite bath and for times shorter than the recurrence time, 
C{t) exhibits the typical damped oscillating behavior of the continuum limit, apart from a shift 
in the oscillation frequency Q. On the other hand, cr'^j^{t) and o'py(i) exhibit in the finite-size 
case the same dissipative behavior as in the continuum case, with a characteristic time l/F. 
However, while (Tprp{t) seems to reach, before the Poincare recurrence time, the same plateau 
value k^T as in the continuum case, cr'^xi^) appears to reach a value lower than the one expected. 




Figure 5. Plot of Cn{uj) vs. uj for different values of the size N of the bath and of the maximal 
frequancy ojc (center) or the width oq, for a uniform (top), Lorentzian (center) and adjusted 
Lorentzian (bottom) sampling of bath frequencies. The continuum limit C{uj) corresponds to 
the dashed line. 



i.e., l/{j3^1). These effects are due to the fact that the frequency shift Yle^l/'^'i is different 
from the continuous limit one na, which appears in J^q. In fact [/3 (ujq — e|/a;|)] is equal 
to the plateau value of cr'^T^t) reached before the recurrence. 

In the over-damping the effect of the frequency shift is larger, one observes a larger 

difference between the continuum and the finite case. In fact, while the continuum limit curves 
display the typical over-damped behavior without any oscillations, the finite-case curves exhibit 
the same behavior observed in the under-damping case. Moreover the difference between the 
plateau values before the recurrence for o"^ j,(i), is also much larger. 

It is clear from Figure |4] that with the Lorenzian sampling one obtains curves that behave 
more similarly to the continuum ones, for the same bath size and recurrence times, with respect 
to the uniform case. This holds both for the oscillation frequency of C{t) and the plateau value 
reached by cr'^j^if) before the recurrence. One may notice the optimal agreement of the curves 
obtained with the adjusted Lorentzian distribution with the continuum ones. 

The same qualitative behavior of the finite size frequency sampling appear in the Fourier 
transform of C{t). In Figure [5| with the same parameters of Figure [4| one notices that in the 
under-damping regime C{uj) is characterized by a peak corresponding to the oscillation frequency 
VL. a similar curve characterizes Xq^^/ {z1,/S.y), but the position of the peak is shifted. This shift 



corresponds to the change in the oscillation frequency of Ci\i{t) with respect to C{t). In the 
over-damping case the = 600 curve maintains the look of the under-damping case, while the 
continuous one looses the peak, then confirming that in this case there is a worse agreement 
between the continuous and the finite cases. 



T = 1000, K. = 0.09, a = 10.0, loq = 1.0 



T = 0.001, K = 0.09, a = 10.0, = 1-0 




Figure 6. Entropy production vs. time t, according to the the three definitions AS*^, AS'^'' 
and ASi, for a uniform distribution of bath frequencies with a cut-off Uc = 30, for different sizes 
N: 3 (blue), 25 (green), 600 (red), and for the adjusted Lorentzian distribution with oq = 0.1 
(yellow) for = 600. The dashed black line corresponds to the continuum limit. The initial 
conditions are those of Figure [2] in the classical case, apart from the first moments (Q^(0)) = 10 
and (P^(0)> = 10, while they are those of Figure p3l in the quantum case. 



Things improve when ujc and the size become larger, keeping the frequency density 
constant. In this case the peak shifts towards its continuum position in the under-damping 
case, while in the over-damping case the peak tends to disappear. This improvement is due 
to the fact that the frequency shift Ylii^'i/'^'i approaches doj '^{vj) /uj^ = Ka. In fact the 
difference between these quantities is due to two terms: one given by the difference between 
the sum X^^ef/wf and the integral up to Wc, which is of order 1/A^ and is negligible for the 
sizes reported in Figure [5} and one, more relevant, corresponding to the contribution to the 
integral arising from frequencies larger than lOq. This term is proportional to k. Thus, in order 
to maintain the difference between X^^^l/"^! ^"^^ constant, ujc must increase as k increases. 
In particular for a given set of parameters, which would correspond to over-damping in the 
continuum limit, one would never obtain over-damping behavior if oo^ is too small. 



If ojc or gq are kept fixed, and increases, the behavior remains the same, only the recurrence 
time tp increases and the smaUest frequency zi decreases. 

With the Lorentzian sampling of parameters convergence improves both in the under-damping 
and in the over-damping cases. In fact, by choosing oq and N so that the recurrence time is 
of the same order as in the uniform case, the value of the frequency shift is closer to na. This 
is due to the fact that the highest frequency is much larger. One has to exercise some care in 
choosing oq neither too large (in order to have long recurrence times) nor too small (in order to 
avoid too sparse a sampling close to the highest frequency). 

We also show in Figures |4] and [5] the effect of adjusting the coupling with the highest-frequency 
oscillator. The behavior of the continuum is optimally matched with the choice = 0.1 and 
N = 600. 



6.3. Finite-size entropy production 

We report in figure [6l the behavior of the entropy production according to the three definitions, 
i.e., the Poised {AS^: eq. the Breuer {AS^': eq. (es))) and the ELB {ASf. eq. (|7l[), 

for different values of N in the uniform case and for N = 600 for the adjusted Lorentzian 
cases. One notices that in the adjusted Lorentzian case one reaches an almost perfect agreement 
with the continuum limit already for N = 600. In the quantum case the expression of the 
entropy production Ai^* obtained with the adjusted Lorentzian binning does not approximate 
perfectly the continuum limit. This is due to the poo r con vergence of the term {Q(t)r](t)) which 
is contained in the averaged interaction energy (Hj) (F.l). The same can be observed for ASf , 
due to the noisy behavior of A{t) and A{t). At finite sizes the Breuer expression AiS^"^ can 
assume negative values, whereas both A^S and AiS^ remain positive. However, in the uniform 
case one obtains a slower convergence with respect to the Lorentzian case, both in the adjusted 
and in the non-adjusted case (not shown). 



7. Bath entropy 

The bath entropy at time t is given by 



Ssit) = -trpB(t) lnpB{t), 



(112) 



where pB{t) = trg p{t) is the reduced bath density matrix. Since the total density matrix is not 
a product state ps <8> Pb at times t > 0, one cannot simply split the total entropy into system 
entropy plus bath entropy. Thus one introduces the correlation entropy Sc- 



Si 



tot 



tr pir) In pit) = S{t) + SBit) + Scit). 



(113) 



We note that —Sc{t) is the mutual information between the central oscillator and the bath [27]. 
Since the total entropy is conserved and the initial correlations vanish, one has S'c(O) = and, 
according to this definition, 

S,it) = -ASit) - ASBit). (114) 



We can easily verify that [13] 



Scit) = -D 



Pit) 



At)llpr{t) 



< 0. 



(115) 



Thus the correlation entropy is always negative or zero. By comparing this last equation with 

(116) 



eq. ( 71 ) , one finds jlB] 

AiS{t) + 5c(t) = -(3 (AHb) - ASBit) = D[pBit)\\p'^] > 0. 



A5s Sc D[pB{tW^] 




t t t 



Figure 7. Entropy change in the bath ASb, correlation entropy Sc and Kulbach-Leibler 
distance between psit) and in the high temperature T = 1000 regime, for different values of 
the coupling k, for a uniform sampling of frequencies with ojc = 30 (blue), and for a adjusted 
Lorentzian one with ao = 0.1 (red). The bath size is = 600 and other relevant parameters 
are set as for figure ([6]). 



One notices that if the approximation of a bath remaining at equilibrium (ideal bath) were 
valid, i.e., pB(t) = p^, the correlation entropy would be equal to minus the entropy production 
Scit) = —AiS{t). In this case the variation of the bath entropy would be equal to the heat flow. 

This 



The method to numerically evaluate the bath entropy is detailed in Appendix G 



calculation relies on the fact that the bath density matrix is Gaussian at each time, and therefore 
is fully characterized by the time-evolving bath covariance matrix (Appendix H). We now turn 
to the discussion of the results. 

As observed before for the central oscillator entropy, with an adjusted Lorentzian sampling of 
frequencies the Poincare recurrence times are much longer, and one can observe a convergence 
toward an asymptotic plateau of the different quantities, as shown for the classical case in figure[7| 
In particular the asymptotic value of ASb does not depend on the coupling constant k, while, in 



agreement with (G.IO), it depends on the initial variances of the central system (see Appendix] 



[H| ). Interestingly, in the over-damping case ASb increases at the beginning, reaching a maximum 
independently of the initial conditions, before decreasing to the equilibrium value, much as minus 
the interaction energy (93) does (not shown). The relaxation time is longer than that of the 
central oscillator. This does not happen in the under-damped regime. 




Once the bath entropy ASb is evaluated, o ne al so gets the correlation entrop y 5c using 

using (116). is negative by definition (115), and 



(114) and then the distance D[pb 



its absolute value grows with the coupling as —ln(uj^ — na), since Ka approaches Wq, similarly 
to the entropy of the system AS" (100). Like the interaction energy term (93), the asymptotic 



value of Sc does not depend on the initial conditions of the central system (data not shown). It 
turns out instead that the asymptotic value of 5c vanishes as k — )• 0. This is confirmed by the 
fact that the coupling-independent asymptotic value of A5b equals minus the central oscillator 
entropy change in the limit of vanishing coupling: 



A5b(oo) = - A5(oo)|^_ 



In 



Wo 



InA(O) 



(117) 



The distance D{pB{t)\\p'^) increases with the coupling as Ka/{ujQ — na), like the negative 
entropy flow (cf. equation (|99|)). This quantity does not vanish for k — t- 0, where it equals the 
entropy production: 



K-5-0 



T 



[A5(oo) - Ac5(oo)] 



(118) 



In 



UJo 



InA(O) -/3 



As a consequence, the bath density matrix operator is always changed and the KuUback- 
Leibler distance from the density operator at canonical equilibrium becomes larger when 
increasing the coupling. This suggests that for our model the ideal bath approximation, namely 
the assumption pb(oo) ~ p^, which would imply A5b — — AgS = I3A{Hb) (116), is not 
valid even in the thermodynamic limit. One observes in figure [8] that in the quantum case the 
asymptotic value of the bath entropy change appears to grow with the coupling, which could 
be an effect of the entanglement or quantum correlations between the bath and the system. 
The dependence on the coupling is apparently weaker than that exhibited by the entropy flow. 
This means that also here the Kullback-Leibler distance between psit) and pg' is relevant and 
strongly increases with the coupling. Due to the quantum contribution in the interaction term, 
its asymptotic value can be orders of magnitude bigger than the one assumed in the high- 
temperature limit, in the same way ELB entropy production does. In the limit of vanishing 
coupling analogous considerations of the classical case can be made, since the correlation entropy 
Sr vanishes. 



8. Conclusions 

In this paper we studied the thermodynamic description of a process of transient relaxation in the 
QBM model where a central harmonic oscillator initially prepared in a Gaussian nonequilibrium 
state is bi-linearly coupled with a bath of harmonic oscillators initially prepared at equilibrium. 

We compared two ways of defining entropy production during the ensuing relaxation process 
of the central oscillator. Both definitions are expressed as the difference between the change 
is the von Neumann entropy of the system minus an heat divided by the temperature of the 
reservoir. The 'ELB' one is based on defining this heat as minus the energy change in the 
bath and thus has a straightforward physical interpretation, while the 'Poised' one (beyond 
non-Markovian transients) defines heat in a less transparent way in term of the change in an 
effective "mean force" Hamiltonian. Both expressions are positive by definition but in a general 
non-Markovian quantum regime they both may exhibit oscillations. However, in the Markovian 
limit, while the 'ELB' may still exhibit oscillations, the 'Poised' one becomes a monotonically 
increasing function of time. The two definitions coincide for vanishing coupling but we have 
shown that for finite coupling the 'ELB' is always larger than the 'Poised' one. Their difference 
contains the expectation value of the interaction Hamiltonian and can thus be made arbitrarily 
large. In the low-temperature limit the contribution due to the quantum corrections in the 
interaction term can make this difference order of magnitudes larger than in the classical case. 
Finally, we showed that in the classical over-damped regime the 'Poised' one converges to the 
entropy production defined in stochastic thermodynamics. 

We numerically studied the exact dynamics of our system for a finite number of oscillators 
in the bath. Using two different samplings of the bath frequencies, a uniform and a Lorentzian 
ones. In both cases the period of the Poincare recurrences increases with growing density of 
bath frequencies but the Lorentzian sampling guarantees a faster convergence to the continuum 
limit curves as a function of A^. 

Finally, we numerically studied the evolution of the von Neumann entropy of the bath which 
results from the relaxation process of the central oscillator. This enabled us to calculate the 
evolution of the system-bath correlation entropy (or minus the mutual information) and the 
Kullback-Leibler divergence between the bath density matrix at time t and its initial thermal 
equilibrium form. We observed that for a given initial condition of the central oscillator, the 
asymptotic value of the bath entropy change does not depend on the coupling in the classical 
limit, while it slightly does in the quantum regime. In the limit of vanishing coupling strength 
the correlation entropy vanishes, what means that the change in the von Neumann entropy of 
the bath becomes equal to minus the change in the central system entropy. We also observed 
that the Kullback-Leibler divergence of the bath density matrix never vanishes, thus indicating 



that the assumption of an ideal bath which always remains at equilibrium is not satisfied. As 
expected, this divergence grows significantly with the coupling, as the ELB expression of the 
entropy production does. 

While our study revealed important features in the QBM model, it also indicates that no 
definite formulation of a consistent thermodynamics of out-of-equilibrium quantum systems in 
presence of non-vanishing coupling with the bath is yet available. 
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Appendix A. Equivalence of the UUersma and Fleming solutions 

The convolution with the noise in (15) corresponds to the sum over the bath operators in ^ 
and the expression for bath operators in ^ correspond exactly to that in ^ . This equivalence 
is recovered thanks to the following equation relating the quantities ^io(t) to the propagator 
A{t), and to the equation relating quantities Aij{t) to Ajo(t), where i and j are bath indices: 

Aio{t) + LofAioit) = -eiA{ty, (A.l) 
Aij{t) + Lo]Aij{t) = -ejAioit). (A.2) 

Given the initial conditions AiQ(t) = Aio{t) = Aij{t) = and Aij{t) = 6ij these equations imply 
that 

Aioit) = -ei drAir)^^ ^; (A.3) 

JO 
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Appendix B. Covariance matrix 

In order to evaluate the covariance matrix of the central oscillator, one evaluates the first 
and second moments for position and momentum Heisenberg operators. From the Heisenberg 
solutions Q, one obtains by averaging over initial conditions (18) (20) (23) the following 
expressions: 



m))=m (Q(0))+A(t) (P(0)); 

{p{t)) = A{t) m))+m (p(o)); 

'ylit) = {Q\t)) - {Q{t)f 

= A{tfag{Of + 2A{t)A{t)Cgp{0) + A{tfap{Of + ^^(i); 

= A{tfag{Of + 2A{t)A{t)Cgp{0) + A{t)%{Of + a^^Ct); 
Cgp{t) = l{{Q{t)-{Q{t)),Pit)-{Pm) 



2 

1 d 

2dt 



(B.l) 
(B.2) 

(B.3) 
(B.4) 

(B.5) 



The thermal parts of the covariance matrix for a finite bath have the form 

N 



1=1 

N 



(B.6) 
(B.7) 



The latter can be also generally written in integral form as shown in (33). By combining 



equations (B.6) and (B.l) one easily obtains also the second moments of the momentmn and 



position operators. 



With the continuous bath with the Ullersma coupling strength (37) and the high cut-off limit 



(41), by inserting (43) in the equations (33), one obtains 



<^lAt) = {l + aUt)} (Q2)^^ + Al^{t) + 2 [AUt)C{t) - aioc{t)C{t)] ; (B.^ 



'eq 



+ 2 



(B.9) 



where aioc(t) = A\oc{t) + 2rAioc{t), and the position equilibrium correlation function (B.IO) is 
given by 



Cit) = ^{{Q{t),QiO)})^^ 

j-QO 

/ dt%oc(t')e'"*'- , . 
Jo ■ 



2 



^ 
1 



Jo 



10 (r - iuy + 02 

When calculated with Aioc(t), it differs from the exact one by corrections of O (1/a): 
C{t) = ^aioc(t) + ^ Im {e-(r+^^)* [^{1 + (F + il7)r^) 

1 cot(7rara)e~"* 



(B.IO) 
(B.ll) 



fB.12) 



-y 



. (B.13) 



In this expression, 'ip{z) = dlnrE(z)/dz is the digamma function, and was defined in 46 The 



last term contains the so called thermal transients, which vanish slowly in the low-temperature 
limit. By discarding terms of O (l/a^^^ in the quantum limit or^ ^ 1 it can be approximated 
for 4 3> 1/q by the series 



Im^ 



-et/Tfi 



(B.14) 



One should remark however that its second time derivative diverges at t = 0, and that other 
terms should be taken in account in order to remove this divergence. We truncate this sum to 50 



terms, what guarantees a good description for t ^ 1/a. One has however to take into account 
the fact that our approximations do not describe well the behavior for t < 1/a. 



at t 



The equilibrium second moments one gets from the equilibrium correlation function (B.IO) 
are given by 



oq 



C(0) 







2 



eq 



TT 



C7(0) = 1/(3 + -hV Re (In a/i^ - + (F + in)Tp)) 

h{Q^ - r^) 



+ 



(B.15) 



(B.16) 



Here the average is carried over the equilibrium density matrix p'"^ (55). This time, in order to 
correctly evaluate (-P^)gq to O (l/ct), one has to consider all the terms contained in Ca,u{t)- One 



should keep in mind that (-P^)gq contains a contribution Ina/v, what explains the necessity of 



introducing a high-frequency cut-off. 



Appendix C. Effect of initial slips 

As we have recalled in section ([3]), slips in the averaged momentum operator and in the 



correlation matrix are produced by the kick-like force term F{t) in the QLE (10), acting during 



an inital time interval of duration ~ 1/a ~ 1/A. The local propagator (43) actually contains 
such slips from t = 0^, so that it is correct apart from corrections of O (1/A) only for t 3> 1/A 
(see eq. |50[), namely when the kick vanishes. In fact, while A{0) contains corrections of the kind 
(1/A^)e^^^^ with respect to Aioc(O), its second time derivative contains a term e~^*, so that it is 



negligible only for t ^ 1/A. With the local propagator (43) we are not going to consider the 



detail of the evolution in the initial time interval of duration 1/a, which is considered to be 
much shorter than the other time scales in the large cut-off limit. 
The initial slips correspond to a fast shift of the initial conditions: 



(P(0))^-2r(Q(0)) + (P(0)) 



(p2(o)> ^ (q2(o)\ - 2r 



(W(CWU(no)>; 



(C.l) 



Let us now discuss the effect of the initial slips (C.l) on the Breuer (68) and ELB (71) 



entropy definitions. They appear as a nonvanishing value for lim(^o+ ^cS{t) = AcS'(0'^) and 
limt^o+ ^S{t) = A5(0+). 

For the Breuer entropy flow one has, by using the expressions of the moments reported 



in (B.l): 



1 



(4r(Q2(o))-2rc,p(o)) 



(C.2) 



which is due to the shift on (|p2(0)). For the entropy change one has 

AS{0+ 



A. 



-) = (Ai,(0) + 1) ln(Ai,(0) + 1) - Ais In A^, 
-(A(0) + 1) ln(A(0) + 1) - A(0) In A(0); 



(C.3) 
(C.4) 



where A(0) is given in (85). Ajg > A(0) > satisfies the Heisenberg principle, and therefore 



A5(0"'") is always positive. 

For the ELB entropy flow one has 



Ae5(0+) = /? [4r (q2(o)) _ 2rCgp(0)) - Ka (Q2(0)>] . 



(C.5) 



Here the sign is determined by the last term, which is generally larger than the first one, due to 
the large value assumed by a. 

Appendix D. Calculation of the Poised entropy production 

In order to evaluate the Poised entropy production Ai^^ in eq.(61), we have to find the poised 



density matrix /Og(i). Once it is known, on can evaluate AjS"^ as 

Ai^P = A5 - Ae^P = A5 - Trs (/9s (0) - ps{t)) Inp^(t). 



(D.l) 



In order to evaluate p%,{t), we rewrite equation V{t)p'^{t) = Pg{t) in the Fourier transform 
space associated with the corresponding Wigner W^{q,p,t). Using eq.(34) we get 



Ws*(*"^(t)fc,t)e-5'=^'^^W'^ = W^ik,t). (D.2) 
Since we only consider initial Gaussian distributions, we have seen that the solution remains 



Gaussian at any time. Therefore, to solve (D.2), we look for solutions of the form 



where CT*{t) is a symmetric 2x2 covariance matrix: 



Cqp{t) cr* (t) 



(D.3) 



(D.4) 



and the vector z*{t) contains the first moments q*{t) and p*{t). By using expression (D.3) in 



equation (D.2) one straightforwardly finds the relation between the covariance matrices and the 



first moments: 



^{t)z*{t) = z*{t). 



(D.5) 
(D.6) 



is equivalent to a system of three equations. One gets 



From equation (D.6) one finds z*{t) = at all times. Equation (D.5) for the covariance matrix 

(D.7) 



cut) 



Dn{t)Di2{t)D2i{t) 

Cqp,T{t)S22{t) + alT{t)S2l{t) + t7^,T(0^23(t) _ 

Dn{t)Di2{t)D2i{t) 



Here 



A^it)Ait)Aity, 



(D. 



Suit) = 533 = 1 - 2A\t) + A\t) - A{t)A{t) 
Si2{t) = -2AA{t){l - A\t) + A{t)A{t)y, 
Sisit) = -A\t){l + A\t) - Ait)Ait)); 
52i(t) = i(t)i(t); 
S22{t) = I - A\t) - A{t)A{t); 
S23it) = Ait)A{t); 

S,,{t) = -A\t)il + A\t)-A{t)A{t)); 
532(t) = -2i(t)i(t)(l - A\t) + Ait)Ait)y, 
Dn{t) = 1 - A\t) + A{t)A{ty 
Di2{t) = 1 - 2i(t) + A\t) - A{t)A{ty, 
D2i{t) = 1 + 2i(t) + A^{t) - A{t)A{t). 

Thus the elements of the matrix cT*{t) are expressed as Hnear combinations of thermal covariance 



elements, whose coefficients are functions of the propagator matrix elements. Expressions (D.7) 



are valid both in the finite-size case and in the thermodynamic limit. In order to get the poised 
covariance matrix at t = one has to evaluate the f — )• 0"^ limit of (D.7). 



Appendix E. Liouvillian operator in the adjoint dynamics 

In the adjoint dynamics, the system and bath Heisenberg operators satisfy the following 
equations: 



Q^{t) = -'-[H,Q,{t)] , P^{t) = -'-[H,P^{t)], 



(E.l) 



where one has a change of sign respect to the usual Heisenberg dynamics Q. By proceeding 
in the same way as in the usual case, we find the following modifications involving the matrix 



QLE satisfied by the central oscillator (10): 

H{t) ^ H^'^it) = ^j^^^^ 

N 



Qi{0) coscjjt H — ^ sina;,;t 



(E.2) 
(E.3) 



This imply that the evolution is given as in equation (15), but with a change in the off-diagonal 



elements of the matrix propagator (16): 



A{t) -Ait) 
-A{t) A{t) 



(E.4) 



while the noise kernel (45) remains unchanged. By considering the above modifications, one 



straightforwardly find from (28) the FP-like equation satisfied by the reduced Wigner for the 



central oscillator for the adjoint dynamics: it remains exactly the same except for a change in 
the signs of the drift term, the harmonic forcing term and the anomalous diffusion coefficient. 
This imply that the limits generally considered in section ([s]) for the usual process do also apply 
for the adjoint, moreover the late time density matrix is the same as in ([s]). 



Appendix F. Interaction energy term 

A way to express the average of Hi, which is useful in the continuum hmit, is the following. By 
plugging the Heisenberg formal solutions for the bath operators as functions of Q{t) (|9| in the 
interaction term, one gets 

{Hj{t)) = ^^Q(t)Qi(t)^ = {Q{t)v(.t)) + (^l^s k{t - s)Q{t)Q{s)^ , (F.l) 

where ri{t) is the fluctuating force term defined in (12). Using ^ and the initial absence of 
system-bath correlation, which implies {Q{0)ri{t)) = 0, one gets 



{QitMt)) = - fdt'A{t-t'){rj{t')v{t)) 
Jo 

^^^^-E{u,T) [ ds A{s) cos{ujs) 
Jo 



do; 



(F.2) 



Interestingly the use of the Ullerma strength (37), from which one obtains that the damping 
coefficients are time-independent in the large cut-off limit, implies exactly the same integral 
form for the anomalous diffusion coefficient, so that {Q(t)r](t)) = Dqp{t). 

In fact in the case of the Ullersma strength (37) with large cut-off the diffusion matrix (30) 
can be approximated by [T7] 



1 



dr 



v{t - r)*"^(t - r) + *(t - T)u'^{t - t) 



(F.3) 



This straightforwardly leads to the equivalence between {Q{t)rj{t)) , which is part of the average 
interaction term (|93[ ), and the anomalous diffusion term Dqp[t). 

The integral ( |F.2[ ) is done by using the local propagator Aioc{t) (43), by first integrating over 

Fcit); 



time and then in the complex u plane. One obtains 



Dgj,{t) = {p^)^^-nl{Q 



' eq 



Alocit) + Alocit) 2T 



d 



dt 



Fc{t) = -{C{t) + nlC{t) + 2TC{t)). 



(F.4) 
(F.5) 



The time dependent term contained in (F.4) vanishes so that one recovers the late-time 
anomalous diffusion coefficient (51), which is a positive quantity. Here we have made use of 
the approximate equilibrium correlation function (B.12). It follows that 

(aT^)2 £e-«/"/3 



Fc{t) = Kh 



1 



cot(7rar^)e' 



-at 



2 (1 + (r/a))2 + (J^/a)2 vr^(ar^)2-£2(^ + rr^)2 + (Or^)2 
which can approximated for t ^ 1/a by 



(F.6) 



Fc{t) 
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' + rr^)2 + (J7r^)2- 



(F.7) 



In the classical limit the anomalous diffusion coefficient vanishes. 

To evaluate the second term of the sum in (F.l), we can use the fact that for large cut-off 
K(t) ~ —Ka5{t). Then one gets for the interaction term complete expression (93). 



In the finite case one has, by using the solutions (|3| and the initial conditions (18): 



N \ N 



Q{t)Qi{t) )=Y.{ ^i^OoMo (Q'(0)) + AooAio (p2(o)) + (ioo^iO + i^o^oo)C,p(0) 



ii=l / i=l 

TV 



This expression is useful for a numerical calculation in the finite case. By using the expression 
of the AiQ{t) and the Aij{t) in function of the propagator A{t), one obtains the expressions 
( F.1|F.2 ) exploited in the continuum limit. 



Appendix G. Evaluation of the bath entropy 

One can straightforwardly evaluate the general quantum bath entropy (112) if one finds a 
coordinate transformation that puts the density operator p-Q in a normal form, namely a product 
of independent oscillator thermal states: 



/>B = (^(l-e-^Oe-^^"^ (G.l) 

e 

where n£ = a\ae, with ae = {qe + ipe)/\^, and where the /?£ are suitable effective inverse 
temperatures. In fact, by putting the density operator in this form, the calculation of entropy is 
easily obtained by carrying the trace over the space of the eigenstates of the number operator: 
|ni, n2, . . . ,n£, . . . , njy). One obtains 

Sb = Y1 (i^e + V2) Hh + 1/2) - ih - 1/2) ln(A;, - 1/2)), (G.2) 

e 

where = ^ coth = ((?|) = {pf) = {ui) + ^. For simplicity we have put here h = 1. 



We know from (2.2) that the reduced density matrix for the bath is Gaussian. Here first 
moments can be shifted to 0, as this transformation leaves the entropy invariant. Then f rom an 
informational point of view the bath is fully characterized by the covariance matrix af^ (H.l). 



The normal form (G.l) and values of the fc^'s can be actually recovered by a "pseudo- 
diagonalization" of the correlation matrix. This can be done using a symplectic transformation, 
I— )• S.^ where /S is a 2N x 2A^-matrix, i.e. a transformation preserving the bosonic commutation 
rules: 

P = S(3S^, (G.3) 

where 

/3=(_°1 J); 1 = (M. A;,^ = l,2,...,iV. (G.4) 
One then choose S such that the correlation matrix in the new basis is diagonal: 

fj^ I—)- cj'^ = Sa^S'^ = diag(Ki, K2, ■ ■ ■ , k-n, k-i, K2, ■ ■ ■ , k^). (G.5) 

This can always be done, as affirmed by Williamson's theorem [28j, due to the fact that the 
correlation matrix is symmetric and positive definite. Due to the particular block form of the 
correlation matrix, the fc^'s are doubly degenerate, as shown in |29| . 

The pseudo-eigenvalues and the symplectic matrix S can be obtained, as explained in [30], by 
diagonalizing the symmetric matrix K f3a^ (3'^ K'^ , where the matrix K is obtained by a Cholesky 
decomposition of the correlation matrix: 



a 



B = k'^K. (G.6) 



This can actually be carried out, since is positive definite. The eigenvalues one finds are 
actually the doubly-degenerate squares of the pseudo-eigenvalues ki. 
Given an operator A, its Wigner transform is defined by [31] 



A 



q + 



(G.7) 



Mp, q)= Jdz e'P'/" (^q - 

In the general quantum case, the Wigner transform of In /3b is —£jM^ — InZ^, since one has 

PB =exp{-e"'MC}/ZB, (G.8) 

where M is a square 2N x 2N matr ix. This matrix transforms under a symplectic transformation 
of the phase space operators (G.3) like (cr^)"^: 



M^M' = {S'^)-^MS-\ 



(G.9) 



In the classical limit the diagonalized matrices M' and (o"'^)^^, coincide and therefore also 
(cr^)^^ and M have to coincide. Therefore in the classical limit the Wigner distribution 
corresponding to the density operator has the same expression as the classical probability 
distribution apart from multiplicative coefficients, i.e.. 



1 



WB(g,p,t) =exp{--et(ai^)-ie}/ [{27Tf{deta'')'/^ 



(G.IO) 



where is given in eq. (H.l). Then the entropy of the bath can be easily calculated via a 
Gaussian integral: 

Sb = - j dqdpWB{q,p,t)ci.\n[{2n)^WB{q,P,t)c^] = N + \u{deia^f'^ . (G.ll) 
This result can be also be obtained by noticing that in the classical limit, where /c£ ^/2, the 



expression (G.2) for the bath entropy reduces to 

N 

Sb = N + \nWkt = N + ln(det Sa^S'^f/^. 



(G.12) 



Thus, since det S = 1 (G.3) we recover (G.ll). 



Appendix H. Bath covariance matrix 

The covariance matrix of the bath is defined as: 



1 



(te,o}>-(c,.)(ei>, 



(H.l) 



where ^ = (Qi, Qn, Pi, Pn), i and j identify the bath oscillators. 

Using the Heisenberg solutions (3]), the variance of a bath position operator with average 
taken over the initial conditions (18) give, with help of (|23|): 



{Qlit)) - {Qi{t)f 



Alaq{Q) + Aio{t)ap{Q) + 2iio(i)^*o(t)Qp(0) 

N 



+E 



Al 



ii 



E{u}e,T). 



(H.2) 



The last sum can be rewritten, using expression (A.4) for the A^is, as foUows: 



N 

E 

e=i 



42 



N 2 ft 
tt^i Jo 



^ E{u}i^T) _^ 2^ /" Aio(r) cos(wir)£'(wj, T). 



(H.3) 



It contains a term expHcitly depending on the initial conditions of the central oscillator, and a 
thermal part. The latter one is made of a term explicitly depending on the initial conditions 
of the bath oscillator, a sum of the kind X]^=i ^^ • ■ - ^ which is easily put into integral form by 
using the strength (14), plus an integral containing Ajo(i) and an oscillating function of time, 
multiplied by the coupling of the oscillator. The same structure is obtained for every term of 
the bath covariance matrix afj between any momentum and bath operators. 
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